home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zbesj.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  4.4 KB  |  123 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((hpi 1.5707963267948966))
  12.   (declare (type double-float hpi))
  13.   (defun zbesj (zr zi fnu kode n cyr cyi nz ierr)
  14.     (declare (type double-float zr zi fnu)
  15.              (type (simple-array double-float (*)) cyr cyi)
  16.              (type f2cl-lib:integer4 kode n nz ierr))
  17.     (prog ((i 0) (inu 0) (inuh 0) (ir 0) (k 0) (k1 0) (k2 0) (nl 0) (aa 0.0)
  18.            (alim 0.0) (arg 0.0) (cii 0.0) (csgni 0.0) (csgnr 0.0) (dig 0.0)
  19.            (elim 0.0) (fnul 0.0) (rl 0.0) (r1m5 0.0) (str 0.0) (tol 0.0)
  20.            (zni 0.0) (znr 0.0) (bb 0.0) (fn 0.0) (az 0.0) (ascle 0.0)
  21.            (rtol 0.0) (atol 0.0) (sti 0.0))
  22.       (declare
  23.        (type double-float sti atol rtol ascle az fn bb znr zni tol str r1m5 rl
  24.         fnul elim dig csgnr csgni cii arg alim aa)
  25.        (type f2cl-lib:integer4 nl k2 k1 k ir inuh inu i))
  26.       (setf ierr 0)
  27.       (setf nz 0)
  28.       (if (< fnu 0.0) (setf ierr 1))
  29.       (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  30.       (if (< n 1) (setf ierr 1))
  31.       (if (/= ierr 0) (go end_label))
  32.       (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  33.       (setf k1 (f2cl-lib:i1mach 15))
  34.       (setf k2 (f2cl-lib:i1mach 16))
  35.       (setf r1m5 (f2cl-lib:d1mach 5))
  36.       (setf k (f2cl-lib:int (min (abs k1) (abs k2))))
  37.       (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  38.       (setf k1 (f2cl-lib:int-sub (f2cl-lib:i1mach 14) 1))
  39.       (setf aa (* r1m5 k1))
  40.       (setf dig (min aa 18.0))
  41.       (setf aa (* aa 2.303))
  42.       (setf alim (+ elim (max (- aa) -41.45)))
  43.       (setf rl (+ (* 1.2 dig) 3.0))
  44.       (setf fnul (+ 10.0 (* 6.0 (- dig 3.0))))
  45.       (setf az (zabs zr zi))
  46.       (setf fn (+ fnu (f2cl-lib:int-sub n 1)))
  47.       (setf aa (/ 0.5 tol))
  48.       (setf bb (* (f2cl-lib:i1mach 9) 0.5))
  49.       (setf aa (min aa bb))
  50.       (if (> az aa) (go label260))
  51.       (if (> fn aa) (go label260))
  52.       (setf aa (f2cl-lib:fsqrt aa))
  53.       (if (> az aa) (setf ierr 3))
  54.       (if (> fn aa) (setf ierr 3))
  55.       (setf cii 1.0)
  56.       (setf inu (f2cl-lib:int fnu))
  57.       (setf inuh (the f2cl-lib:integer4 (truncate inu 2)))
  58.       (setf ir (f2cl-lib:int-sub inu (f2cl-lib:int-mul 2 inuh)))
  59.       (setf arg (* (- fnu (f2cl-lib:int-sub inu ir)) hpi))
  60.       (setf csgnr (cos arg))
  61.       (setf csgni (sin arg))
  62.       (if (= (mod inuh 2) 0) (go label40))
  63.       (setf csgnr (- csgnr))
  64.       (setf csgni (- csgni))
  65.      label40
  66.       (setf znr zi)
  67.       (setf zni (- zr))
  68.       (if (>= zi 0.0) (go label50))
  69.       (setf znr (- znr))
  70.       (setf zni (- zni))
  71.       (setf csgni (- csgni))
  72.       (setf cii (- cii))
  73.      label50
  74.       (multiple-value-bind
  75.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  76.            var-11 var-12)
  77.           (zbinu znr zni fnu kode n cyr cyi nz rl fnul tol elim alim)
  78.         (declare
  79.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9 var-10
  80.           var-11 var-12))
  81.         (setf nz var-7))
  82.       (if (< nz 0) (go label130))
  83.       (setf nl (f2cl-lib:int-sub n nz))
  84.       (if (= nl 0) (go end_label))
  85.       (setf rtol (/ 1.0 tol))
  86.       (setf ascle (* (f2cl-lib:d1mach 1) rtol 1000.0))
  87.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  88.                     ((> i nl) nil)
  89.         (tagbody
  90.           (setf aa (f2cl-lib:fref cyr (i) ((1 n))))
  91.           (setf bb (f2cl-lib:fref cyi (i) ((1 n))))
  92.           (setf atol 1.0)
  93.           (if (> (max (abs aa) (abs bb)) ascle) (go label55))
  94.           (setf aa (* aa rtol))
  95.           (setf bb (* bb rtol))
  96.           (setf atol tol)
  97.          label55
  98.           (setf str (- (* aa csgnr) (* bb csgni)))
  99.           (setf sti (+ (* aa csgni) (* bb csgnr)))
  100.           (f2cl-lib:fset (f2cl-lib:fref cyr (i) ((1 n))) (* str atol))
  101.           (f2cl-lib:fset (f2cl-lib:fref cyi (i) ((1 n))) (* sti atol))
  102.           (setf str (* (- csgni) cii))
  103.           (setf csgni (* csgnr cii))
  104.           (setf csgnr str)
  105.          label60))
  106.       (go end_label)
  107.      label130
  108.       (if (= nz -2) (go label140))
  109.       (setf nz 0)
  110.       (setf ierr 2)
  111.       (go end_label)
  112.      label140
  113.       (setf nz 0)
  114.       (setf ierr 5)
  115.       (go end_label)
  116.      label260
  117.       (setf nz 0)
  118.       (setf ierr 4)
  119.       (go end_label)
  120.      end_label
  121.       (return (values nil nil nil nil nil nil nil nz ierr)))))
  122.  
  123.